% fields H for rho > a & rho < a

function [Hx, Hphi, Hz, Erho, Ephi, Ez] = func_HE_inside(y, nx, CONSTS, coef, p1)

    k0 = CONSTS.k0;
    eta   = CONSTS.eta;
    eps   = CONSTS.eps;
    g     = CONSTS.g;
    
    B1 = coef(1);
    B2 = coef(2);

    [nxy1 nxy2] = func_nxyk_from_nz(nz, CONSTS);
    [ny1 ny2] = func_nyk_from_nxyk(nz, nx, CONSTS);
    [n1 n2] = func_nk_from_nz(nz, CONSTS);
    [alfa_1 alfa_2 beta_1 beta_2] = func_alfa_beta_from_nz(nz, CONSTS);

    Hx = -n1.*B1.*((-beta_1.*nx + 1i*ny1)./(nxy1)).*exp(-1i*k0*ny1.*y) + ...
         n2.*B2.*((beta_2.*nx + 1i*ny2)./(nxy2)).*exp(-1i*k0*ny2.*y);
     
    Hy = -1i*
    
    Hrho  = -1i*(B1*(p1*bessel_mp1_q1 - n1*beta_1*m*(bessel_m_q1./(Q1*r))) + ...
                         B2*(p1*bessel_mp1_q2 - n2*beta_2*m*(bessel_m_q2./(Q2*r))));
    
    Hphi  = -(B1*n1*(bessel_mp1_q1 - beta_1*m*(bessel_m_q1./(Q1*r))) + ...
                      B2*n2*(bessel_mp1_q2 - beta_2*m*(bessel_m_q2./(Q2*r))));
    
    Hz    = -(B1*q1*bessel_m_q1 + B2*q2*bessel_m_q2);
    

    Erho  = -(B1*(((n1*p1+g)/eps)*bessel_mp1_q1 + alfa_1*m*(bessel_m_q1./(Q1*r))) + ...
                      B2*(((n2*p1+g)/eps)*bessel_mp1_q2 + alfa_2*m*(bessel_m_q2./(Q2*r))));
    
    Ephi  = 1i*(B1*(bessel_mp1_q1 + alfa_1*m*(bessel_m_q1./(Q1*r))) + ...
                        B2*(bessel_mp1_q2 + alfa_2*m*(bessel_m_q2./(Q2*r))));
    
    Ez    = (1i/eta)*(B1*n1*q1*bessel_m_q1 + B2*n2*q2*bessel_m_q2);
    
end
